数学模型

Wang Haihua

🍈 🍉🍊 🍋 🍌

GM(2,1)模型

GM(1,1)模型适用于具有较强指数规律的序列，只能描述单调的变化过程，对于非单调的摆动发展序列或有饱和的S形序列，可以考虑建立GM(2,1)，DGM和Verhulst模型。下面我们介绍介绍GM(2,1)模型。

定义1 设原始序列 $$ x^{(0)}=\left(x^{(0)}(1), x^{(0)}(2), \cdots, x^{(0)}(n)\right), $$ 其 1 次累加生成序列 (1-AGO) $x^{(1)}$ 和 1 次累减生成序列 (1-IAGO) $\alpha^{(1)} x^{(0)}$ 分别为 $x^{(1)}=\left(x^{(1)}(1), x^{(1)}(2), \cdots, x^{(1)}(n)\right)$, 和 $\quad \alpha^{(1)} x^{(0)}=\left(\alpha^{(1)} x^{(0)}(2), \cdots, \alpha^{(1)} x^{(0)}(n)\right)$, 其中 $\quad \alpha^{(1)} x^{(0)}(k)=x^{(0)}(k)-x^{(0)}(k-1), k=2,3, \cdots, n$. $x^{(1)}$ 的均值生成序列为 $z^{(1)}=\left(z^{(1)}(2), z^{(1)}(3), \cdots, z^{(1)}(n)\right)$, 则称 $$ \alpha^{(1)} x^{(0)}(k)+a_{1} x^{(0)}(k)+a_{2} z^{(1)}(k)=b $$ 为 $\operatorname{GM}(2,1)$ 模型。

定义 2 称 $$ \frac{d^{2} x^{(1)}(t)}{d t^{2}}+a_{1} \frac{d x^{(1)}(t)}{d t}+a_{2} x^{(1)}(t)=b $$ 为 $G M(2,1)$ 模型的白化方程。

定理 1 设 $x^{(0)}, x^{(1)}, \alpha^{(1)} x^{(0)}$ 如定义 $15.2$ 所述, 且 $$ B=\left[\begin{array}{ccc} -x^{(0)}(2) & -z^{(1)}(2) & 1 \\ -x^{(0)}(3) & -z^{(1)}(3) & 1 \\ \vdots & \vdots & \vdots \\ -x^{(0)}(n) & -z^{(1)}(n) & 1 \end{array}\right], \quad Y=\left[\begin{array}{c} \alpha^{(1)} x^{(0)}(2) \\ \alpha^{(1)} x^{(0)}(3) \\ \vdots \\ \alpha^{(1)} x^{(0)}(n) \end{array}\right]=\left[\begin{array}{c} x^{(0)}(2)-x^{(0)}(1) \\ x^{(0)}(3)-x^{(0)}(2) \\ \vdots \\ x^{(0)}(n)-x^{(0)}(n-1) \end{array}\right], $$ 则 $G M(2,1)$ 模型参数序列 $\boldsymbol{u}=\left[a_{1}, a_{2}, b\right]^{T}$ 的最小二乘估计为 $$ \hat{\boldsymbol{u}}=\left(\boldsymbol{B}^{T} \boldsymbol{B}\right)^{-1} \boldsymbol{B}^{T} \boldsymbol{Y} . $$

案例

问题

已知 $x^{(0)}=(41,49,61,78,96,104)$, 试建立 GM(2,1)模型。

建模

$x^{(0)}$ 的 1-AGO 序列 $x^{(1)}$ 和 1-IAGO 序列 $\alpha^{(1)} x^{(0)}$ 分别为 $$ x^{(1)}=(41,90,151,229,325,429), \alpha^{(1)} x^{(0)}=(8,12,17,18,8), $$ $x^{(1)}$ 的均值生成序列 $z^{(1)}=(65.5,120.5,190,277,377)$, 记 $$ B=\left[\begin{array}{ccc} -x^{(0)}(2) & -z^{(1)}(2) & 1 \\ -x^{(0)}(3) & -z^{(1)}(3) & 1 \\ \vdots & \vdots & \vdots \\ -x^{(0)}(6) & -z^{(1)}(6) & 1 \end{array}\right]=\left[\begin{array}{ccc} -49 & -65.5 & 1 \\ -61 & -120.5 & 1 \\ -78 & -190 & 1 \\ -96 & -277 & 1 \\ -104 & -377 & 1 \end{array}\right], Y=\left[\begin{array}{c} 8 \\ 12 \\ 17 \\ 18 \\ 8 \end{array}\right] $$

则有 $$ \hat{\boldsymbol{u}}=\left[\begin{array}{c} \hat{\boldsymbol{a}}_{1} \\ \hat{\boldsymbol{a}}_{2} \\ \hat{\boldsymbol{b}} \end{array}\right]=\left(B^{T} B\right)^{-1} B^{T} Y=\left[\begin{array}{c} -1.0922 \\ 0.1959 \\ -31.7983 \end{array}\right], $$ 故得 $\operatorname{GM}(2,1)$ 白化模型 $$ \frac{d^{2} x^{(1)}}{d t^{2}}-1.0922 \frac{d x^{(1)}}{d t}+0.1959 x^{(1)}=-31.7983 . $$

利用边值条件 $x^{(1)}(1)=41, x^{(1)}(6)=429$, 解之得时间响应式为 $$ x^{(1)}(t)=203.85 e^{0.22622 t}-0.5325 e^{0.86597 t}-162.317, $$ 于是 $\operatorname{GM}(2,1)$ 时间响应式 $$ \hat{\boldsymbol{x}}^{(1)}(\boldsymbol{k}+\mathbf{1})=\mathbf{2 0 3 . 8 5} \boldsymbol{e}^{0.22622 k}-0.5325 e^{0.86597 k}-162.317 $$ 所以 $$ \hat{\boldsymbol{x}}^{(1)}=(41,92.0148,155.1561,232.3672,324.5220,429) \text {. } $$ 做 IAGO 还原, 有 $$ \hat{\boldsymbol{x}}^{(0)}=(41,51.0148,63.1412,77.2111,92.1548,104.4780) . $$

计算结果见下表。

序号 实际数据 x^((0)) 预测数据 hat(x)^((0)) 残差 x^((0))- hat(x)^((0)) 相对误差 delta(k)
2 49 51.0148 -2.0148 4.1%
3 61 63.1412 -2.1412 3.5%
4 78 77.2111 0.7889 1.0%
5 96 92.1548 3.8452 4.0%
6 104 104.4780 -0.4780 0.5%

代码

import numpy as np
from sympy import Function, diff, dsolve, symbols, solve,exp
x0=np.array([41, 49, 61, 78, 96, 104])
n=len(x0); x1=np.cumsum(x0)  #计算1次累加序列
ax0=np.diff(x0)  #计算一次累减序列
B=np.c_[-x0[1:],-z,np.ones((n-1,1))]
u=np.linalg.pinv(B).dot(ax0)
p=np.r_[1,u[:-1]]  #构造特征多项式
z=0.5*(x1[1:]+x1[:-1])  #计算均值生成序列
r=np.roots(p)  #求特征根
xts=u[2]/u[1]  #常微分方程的特解
c1,c2,t=symbols('c1,c2,t'); eq1=c1+c2+xts-41;
eq2=c1*np.exp(5*r[0])+c2*np.exp(5*r[1])+xts-429
c=solve([eq1,eq2],[c1,c2])
s=c[c1]*exp(r[0]*t)+c[c2]*exp(r[1]*t)+xts  #微分方程的符号解
xt1=[]
for i in range(6): xt1.append(s.subs({t:i}))
xh0=np.r_[xt1[0],np.diff(xt1)]
cha=x0-xh0  #计算残差

delta=np.abs(cha)/x0  #计算相对误差
print(xt1,'\n------------\n',xh0,'\n------------\n',
      cha,'\n--------------\n',delta)

符号解法：

import numpy as np
import sympy as sp
x0=np.array([41,49,61,78,96,104])
n=len(x0)
lamda=x0[:-1]/x0[1:]  #计算级比
rang=[lamda.min(), lamda.max()]  #计算级比的范围
theta=[np.exp(-2/(n+1)),np.exp(2/(n+1))] #计算级比容许范围

x1=np.cumsum(x0)  #累加运算
z=0.5*(x1[1:]+x1[:-1])
B=np.vstack([-x0[1:],-z,np.ones(n-1)]).T
u=np.linalg.pinv(B)@np.diff(x0)  #最小二乘法拟合参数
print("参数u：",u)
sp.var('t'); sp.var('x',cls=sp.Function)  #定义符号变量和函数
eq=x(t).diff(t,2)+u[0]*x(t).diff(t)+u[1]*x(t)-u[2]
s=sp.dsolve(eq,ics={x(0):x0[0],x(5):x1[-1]})  #求微分方程符号解
xt=s.args[1]  #提取解的符号表达式

print('xt=',xt)
fxt=sp.lambdify(t,xt,'numpy')  #转换为匿名函数
yuce1=fxt(np.arange(n))  #求预测值
yuce=np.hstack([x0[0],np.diff(yuce1)])  #还原数据
epsilon=x0-yuce[:n]  #计算已知数据预测的残差
delta=abs(epsilon/x0)  #计算相对误差
print('相对误差：',np.round(delta*100,2))  #显示相对误差